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ABSTRACT 

For  nearly  two  decades  we  have  witnessed  an  intensive  development 
of  a statistical  methodology  for  assessing  length  of  life  and  relia- 
bility of  performance  from  empirical  data.  The  initial  stimulus  for 
research  on  statistical  problems  in  life  testing  and  reliability  Ccuie 
from  the  need  to  answer  pressing  practical  questions  which  could  not  be 
treated  by  the  existing  statistical  techniques.  Because  life  and  per- 
formance tests  are  so  time  consuiiiing  and  expensive  to  run,  it  is  a 
practical  necessity  to  terminate  them  as  soon  as  possible. 

For  tlie  statistician  this  means  developing  estimation  and  decision 
procedure  for  data,  which  are  severely  curtailed  in  one  way  or  another 
long  before  all  items  on  test  have  actually  failed.  The 
estimation  is  more  con^licated  when  the  data  are  truncated,  i.e.  when 
the  observer  loses  track  of  some  individuals  before  death  occur.  The 
product  limit  method  of  Kaplan  and  Meier  is  one  way  of  estimating  p{t) 
when  the  mechanism  causing  truncation  is  independent  of  the  mechanism 
causing  death. 

This  paper  proposes  alternative  estimators  and  compares  them  to 
the  product  limit  method.  A computer  simulation  is  used  to  generate 
the  times  of  death  and  truncation  from  a variety  of  assumed  distribu- 
tions. No  single  estimator  gives  the  best  fit  to  the  "true"  distribu- 
tion of  death  under  all  situations.  However,  other  estimators  are 
shown  to  be  better  than  the  product  limit  estimator  under  all  of  the 
assumed  situations. 
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I.  INTRODUCTION 


Let  the  random  variable  T denote  the  time  that  elapses  until  an 
event  occurs;  the  event  may  for  example  be  an  equipment  failure,  an 
individual's  death,  or  the  detection  of  a target.  Denote  by  p(t)  the 
probability  of  survival  to  time  t, 

p{T  > t}  = p(t) 

Picturesquely,  T is  called  a lifetime,  and  p(t)  is  a survival  pro- 
I bability; 

F(t)  = 1 - p(t)  is  the  distribution  function  of  T. 

In  the  medical  field,  one  might  wish  to  estimate  the  probadsility , 
p(t)  that  a patient  survives  t after  a certain  surgical  procedure 
for  cancer.  In  electronics,  one  wishes  to  estimate  the  probability  of 
continuous  failure-free  operation  of  an  equipment  for  time  t.  In  the 
military,  one  might  be  interested  in  the  probability  of  conducting  a 
certain  mission,  under  specified  environmental  conditions,  without  de- 
tection by  the  enemy.  The  event  of  interest  may  be  a human  death, 
equipment  malfunctions,  or  sonar  detection.  Following  Kaplan  and  Meier, 
Reference  (1) , this  paper  will  refer  to  the  event  of  interest  as  a 
"death".  The  test  element  in  the  sample  may  be  a human,  a radio,  or 
a submarine.  This  paper  will  refer  to  the  test  elements  as  "individuals". 
Suppose  that  observed  values  of  T are  t^ , t , t ,...t  , so  that  N 
lifetimes  are  observed.  In  this  case  an  appropriate  (unbiased)  esti- 
mates of  survival  to  time  t is 
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number  of  t. 's  > t 

p(t)  = i 

w 

Under  many  circumstances  complete  lifetimes  are  not  observed;  censoring 
occurs  at  certains,  , beyond  which  the  life  of  an  individual  is  not 
known.  In  such  cases  construction  of  an  appropriate  estimate  of  the 
survival  probability  is  more  difficult.  In  this  paper  various  estimates 
of  survival  probability  are  studied  when  lifetimes  are  randomly  censored. 
This  means  that  censoring  times  are  assumed  to  be  realizations  random 
variables  independent  of  the  actual  lifetimes. 

The  product-limit  estimator  of  Kaplan  and  Meier,  Reference  (1),  is 
an  accepted  method  of  dealing  with  the  problem  of  censored  data.  This 
paper  presents  thirteen  non-parametric  estimators,  including  the  product 
limit  function.  Censored  data  sets  are  simulated.  The  thirteen  esti- 
mators are  compared  by  examining  their  performance  on  the  simulated  data 
bases. 
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II.  THEORY 


Ihere  are  two  approaches  to  the  empirical  estimation  of  the  survival 
probability,  p(t): 

(1)  one  may  use  the  observed  fraction  of  survivors  at  arbitrarily 
selected  times  (step  function  estimator) , or 

(2)  one  may  focus  attention  on  the  times  of  the  observed  deaths 
(point  estimator). 

The  initial  discussion  is  based  on  the  assumption  that  all  observa- 
tions are  complete,  i.e. , it  is  assumed  that  all  individuals  remain  under 
observation  until  their  time  of  death.  This  initial  assumption  is  for 
the  purpose  of  simplifying  the  discussion.  Then,  later  in  this  paper, 
the  discussion  is  broadened  to  include  incomplete  data  with  observations 
of  both  death  and  censoring  events. 

Survival  Probabilities;  No  Censoring 

~ ^0  ^ ^1  ^ ^2  ' " ‘ ^i  ^'i+l  ^ ‘ ^ sequence  of  fixed 

times.  Then  if  T is  a lifetime 

p(t^)  = p{T  > t^} 

and  denote  the  conditional  probability  of  survival  to  time  t. , given 

survival  to  t,  , by 
i+1  ^ 

P(ti|ti_i)  - p{T  > t^|T  > t^_^} 
p{T  > t^}  p(t^) 

’ pIt  > t._^}  ° p(t^_^) 
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If  ^)  = 0 , define  = 0 . 

Then 


p(t.)  = p(tjt._^)p(t._^)  = P(tJt._^)-p{t._Jt._2)p(t._2) 


= Tr  p(tjit._^) 

j=i 


(2) 


where  p(t^|t^)  = p(t^)  ; p(o)  = 1 . 


Observations  on  Uncensored  Data  at  Fixed  Times 


Let  a sample  of  N individuals  come  under  observation.  They  are 
all  observed  from  birth  (or  the  appropriate  event  defining  tirse  zero) 
until  death.  v;ith  the  first  approach,  preselects  a series  of  times. 


0 < tj^  < t^  < ...  before  examining  the  observed  time  of  death.  In  the 


medical  follow-up  example,  one  might  select  the  times  corresponding  to 
exactly  1,2,3,...  years  after  a surgical  procedure  for  cancer.  An  esti- 


mate of  the  conditional  probability  of  survival  to  t^,  given  survival 


to  t.  , is 
1-1 


P(tjt._^) 


tJ.  - r. 
1 1 


N. 

1 


(3) 


With  elements  were  present  at  the  beginning  of  the  interval,  i.e.. 


at  time  and  r^  elements  failed  during  the  interval. 


For  a set  of  data  which  is  not  censored,  N.  = N.  , - r.  , . Now 

1 1-1  1-1 
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N-r  N-r  -r 
(^)  ( ..  _ -)  - 


N N-r, 


N-r  -...-r.  N-r,-...-r. 

, 1 ,1  1 . 
•“  z I irr— ) 


N-r-...-r.  'N-r  . .-r . , 
1 1-2  1 1-1 


= 1 - 


izi 


Now  the  estimate  p(t^)  is  of  the  form 


P(t^)  = 


N-(r^  + r^  + ...  + r^) 


and  this  is  the  Scime  as 


P(ti)  = ^ 


where  is  the  number  of  the  original  group,  of  size  N,  that  survive 


to  t. . If  it  is  assumed  that  the  N individuals  each  have  the  survival 
1 


probability  p(t),  and  that  they  die  independently,  then  S^,  the  random 


number  that  survive  to  tine  t^  is  binomially  distributed,  with  being 


a realized  value  of  S^.  Then,  considering  the  estimate  as  a random 


variable. 


P(ti)  = ^ 


and 


N p(t  ) 

E[p(ti)]  = — = p(t.) 


and 


Var[p(t^)]  = 


p(t.) (l-p(t.) 


Consequently  p(t^)  is  an  unbiased  and  consistent  estimate  of  p(t^) 


This  is  true  for  every  t^,  and  can  be  shown  to  be  true  for  all  t^,  i=l. 


2, . . .1,  as  will. 
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All  of  this  indicates  that  the  estimate  suggested  is  likely  to  be  a 
good  one  if  the  sample  size,  N,  is  large. 

Clearly  p{t^)  • The  survival  probability,  p{t),  is  thus 

estimated  at  a fixed  sequency  of  times.  At  each  time  point,  t^  being  a 
typical  one,  there  are  r^  fewer  survivors  than  at  t^  where  r^  = 0,1, 
2,...,R.  Consequently  a plot  of  p(t^)  shows  a non-decreasing  step 
function,  with  downward  steps  of  varying  sizes  at  t^jt^,...  . 

If  the  cibove  times  are  close  together,  and  if  the  time  of  death  T, 
has  a density  function,  then  one  can  anticipate  seeing  values  of  r^  that 
are  either  zero  or  unity. 

The  so-called  second  approach  is  really  a limiting  case  of  the  first, 
as  the  time  of  intervals  of  measurement  decrease  indefinitely.  Thus  when 
a death  (or  loss)  occurs  it  is  only  a single  event. 

When  no  losses  take  place,  the  case  now  considered,  the  time  t^  of 
the  ith  death  is  a really  a realization  of  a random  variable,  denoted  by 
t^;  this  means  that  p(_t^)  the  probability  of  surviving  is  a random 
variable.  It  can  be  shown  that  the  expected  value  of  p(^j^)  is 

E[p(t^)]  = , 1=1,2, ...,N 


where  t,  < t < < t . 

— i — Z — N 

The  derivation  involves  integrating 


Elp(t^)] 


00 

/ 


p(t) ' 


N! 


(i-1) ! (N-i+1) 


[l-p(t)]^"^(- 


dp(t) 

dt 


•)  [p(t)] 


N-i 


N-i+l 

N+1 


by  transformation  from  p(t)  to  x;  see  Cramer,  Mathematical  Methods  of 
Statistics,  h.  Cramer,  Princeton  University  Press,  1946. 
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Thus  one  is  led  to  use 


I 

i 


p(t^) 


N-i+l 

N+1 


(4) 


as  an  estimate  of  the  value  of  p(t^),  t^  being  the  ith  time  of  death. 
Expression  (4)  provides  estimator  of  the  survival  function  at  times  of 
observed  deaths  when  there  are  no  losses  because  of  censoring.  The 

estimator  at  the  points  t^:  ^ ^2  ^ ^N'  connected  by 

straight  lines,  or  a step  function  with  step  sizes  1/(N+1)  may  be  used. 

The  estimators  of  equation  (4)  give  intuitively  acceptable  results. 
For  example,  if  the  sample  consists  of  only  a single  individual  (N=l) , 
then  death  is  equally  likely  to  occur  before  or  after  the  time  at  which 
the  true  (but  unknov/n)  survival  function  equals  one  half.  Thus,  the 
result  of  equation  (4)  is  reasonable: 

E[p{t^)]  = I 

The  point  estimates  of  the  second  approach  always  occur  at  the  times  of 
discontinuity  forestimates  from  the  first  approach.  For  example,  con- 
sider a data  base  (H=4)  with  deaths  observed  at  times  1,3,4  and  7.  The 
first  approach  gives  the  following  step  function  estimate  of  the  survival 

1.0  0 £ t < 1 

0.75  1 < t < 3 

0.5  3 < t < 4 

0.25  4 < t < 7 

0.0  t < 7 
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The  second  approach  gives  the  following  point  estimates 


p(o)  = 1.0 
p(l)  = 0.8 
p(3)  = 0.6 
p(4)  = 0.4 
p(7)  = 0.2 

A graphic  comparison  of  the  results  of  the  two  approaches  is  given 
below: 


It  is  difficult  to  decide  how  to  smooth  out  the  step  functions 
that  result  from  the  first  approach.  By  connecting  the  tops  of  the 
"stairsteps, "one  places  an  upper  bound  on  reasonable  estimates.  By 
connecting  the  bottom  corners  of  the  stairsteps,  one  places  a lower 
bound  on  reasonable  estimates.  One  might  draw  a smooth,  decreasing 
curve  that  passes  through  all  (or  almost  all)  of  the  vertical  faces 
of  the  step-function  estimate.  The  second  approach  suggests  method 
of  selecting  a unique  point  on  each  of  these  vertical  segments. 


13 


Incomplete  observations 

When  some  of  the  observations  are  incomplete,  equation 
(4)  requires  modification.  The  expected  value  of  the  survival  function 
at  the  time  of  the  first  observed  death  may  be  written: 

(5) 

Here  is  tlie  effective  size  of  the  sample  during  the  interval  termi- 
nated by  the  time  observed  for  the  first  death  (o,t^) . In  the  special 
case  of  no  censoring  events,  the  value  of  is  unambiguous.  It  is 
equal  to  the  initial  scimple  size  = N)  . In  this  case  equation  (5) 
reduces  to  equation  (4). 

Svibsequent  point  estimates  for  t^,  t^,...  may  be  calculated 
iteratively: 

1 

where  t^  = 0 and  is  the  effective  sample  size  over  the  time  interval 
(t^_^,  t^) . Th"s, 

= IT  (6) 

J=1  1 

Variance  of  the  estimators 

Kaplan  and  Meier,  reference  (1),  give  an  expression  for  the  exact 

calculation  of  the  variance  of  step  functions.  They  also  discuss 

"Greenwood's  formula,"  a large  sample  approximation  that  ignores  terms 
2 

of  order  1/N^  . 
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an  expression  for 


i^eference  (2) , presence  without  derivation 
the  variance  of  estimates  using  the  second  approach  (point  estimators) : 

V(t.)  = var  {E[p(t.)]}  = TT  - TT 

3=1  : j=i 

The  notation  here  follows  that  for  the  estimating  equation  (6) . 
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III.  THE  ESTIMATORS 

This  section  describes  the  nine  non-parametric  estimators  and  four 
jackknife  estimators  of  the  survival  probability.  It  also  describes 
the  parametric  estimator  for  an  exponential  decay  function.  Exponential 
life  distributions  are  the  starting  point  for  much  of  reliability  theory 
and  practice.  The  estimator  derived  from  the  exponential  is  regarded  as 
"par"  when  the  simulated  data  is  based  on  an  underlying  exponential  decay 
distribution  for  deaths.  Thus,  when  deaths  are  exponentially  distributed, 
the  non-parametric  estimators  nay  be  compared  relative  to  each  other, 
and  they  may  be  ccxnpared  with  the  parametric  estimator  as  a standard. 

A hypothetical  data  base,  consisting  of  five  individuals,  is  used 
to  illustrate  each  of  the  estimators.  This  sample  data  base  is  as 
follows: 


Individual 

Time  of  Death 

Time  of  Truncation 

A 

1 

- 

B 

Unknown  (>2) 

2 

C 

3 

- 

D 

Unknown  (>6) 

6 

E 

7 

.. 

The  data  have  been  arranged  in  time  sequence  of  the  death  and  trunca- 
tion events.  In  the  medical  example,  the  data  might  indicate  that 
patients  A,  C and  E were  observed  to  die  exactly  1,  3 and  7 years,  re- 
spectively, after  their  surgery'.  However,  B and  D moved  away  or  other- 
wise became  unavailable  to  the  observer  at  these  times.  Further,  the 
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cause  of  the  unobservability  is  unrelated  to  the  patient's  health  and 
life  expectancy. 

A.  STEP-FUNCTION  ESTIMATORS 

1.  The  First  Estimator,  "p^(t)" 

p^(t)  is  a naive  estimator;  it  is  expected  to  perform  poorly 
relative  to  the  other  estimators,  p^  only  depends  on  the  data  from  in- 
dividuals whose  deaths  are  observed.  It  ignores  any  information  from 
the  partial  lifetimes  noted  for  the  censored  observations.  Pj^(t)  is 
simply  the  fraction  of  individuals  surviving  to  at  least  time  t among 
those  individuals  whose  time  of  death  is  known.  It  is  a step  function: 

t Pl^^> 

0-1  1.0 

1-3  0.667 

3-7  0.333 

7-00  0.00 

The  naive  estimator,  p^(t),  takes  no  account  of  the  successful 
survival  intervals  observed  for  the  censored  individuals.  Therefore 
it  is  biased  in  a downward  (pessimistic)  direction. 

2.  The  Second  Estimator,  "p2(t)" 

P2(t)  is  the  product-limit  estimate.  Kaplan  and  Meier,  refer- 
ence (1) , have  shown  that  this  is  the  maximum  likelihood  estimator.  The 
observed  events,  both  deaths  and  truncations,  are  arranged  in  increasing 
order  of  occurrence;  t^,  t2,...,tj^;  where  H is  the  number  of  individuals 
in  the  scimple. 

Let  p(t^)  denote  the  cumulative  probability  of  survival  of  an 
individual  from  time  zero  to  time  t^.  Let  p(t|t^)  denote  the  conditional 
probability  of  surviving  to  time  t(>  t^) , given  that  the  individual  has 
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already  survived  to  time  t^.  Then, 


P2(ti)  = P2(t._^)  • P2(tjt._^) 


If  we  define  t =0  and  p(0)  =1,  then 
0 


^2^’^i^  - lTP2(tj|t  ) 

:=i 


The  product  limit  estimator  is  in  the  form  of  equation  (E-2)  with 


‘iL.i 


If  the  event  at  t . is 
truncation  ^ 


If  the  event  at  t.  is 
a death  ^ 


Here  n^  is  the  number  of  individuals  observed  surviving  in  the  interval 
tj_^  < t < t j . This  formulation  causes  the  product  limit  estimator  to 
be  insensitive  to  the  exact  time  of  the  censoring  events. 

The  estimator  is  unity  from  time  zero  to  the  time  of  the  first 
event,  t^^,  reflecting  the  fact  tiiat  all  individuals  in  our  example  are 
observed  to  live  until  at  least  time  t^. 

- If  the  event  at  time  t^  is  a truncation,  then  the  estimator 
remains  at  unity  until  at  least  time  t^.  Again,  no  deaths 
are  observed  in  the  sample  before  t^. 

- If  the  event  at  time  t^  is  a death,  then  the  estimator  drops 
to  (N-D/tj,  This  drop  reflects  the  observed  death  of  1/N  of 
the  survival  sample  just  prior  to  t^. 

Values  of  the  estimator  p^  are  calculated  iteratively  at  successive 
values  of  t^ (i=l,2, . . . ,N) . 
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The  size  of  the  survival  sample  declines  as  truncations  and 
deaths  remove  individuals  from  observation.  For  the  hypothetical  data 
base  listed  above,  one  obtains; 
t 


0-1 

5/5 

S 

1.0 

1-2 

4/5 

= 

0.8 

2-3 

(4/5) 

X 

(3/3) 

= 

0.8 

3-6 

(4/5) 

X 

(2/3) 

= 

0.533 

6-7 

(8/15) 

X 

(1/1) 

= 

0.533 

7-00 

(8/15) 

X 

(0/1) 

2= 

0.0 

The  product-limit  estimator  explicitly  accounts  for  the  sur- 
vival of  these  individuals  (up  to  the  time  of  the  last  death  before 
each  censoriny  event).  Thus,  P2(t)  is  a step  function  with  a value 
that  is  not  less  than  Pj^(t)  for  any  value  of  t.  If  the  sample  contains 
no  censoring,  then  Pj^(t)  and  P2(t)  are  identical. 

If  the  last  event  in  the  saitple  is  a truncation  rather  than  a 
death;  then  tlie  modified  data  give  the  following  estimate,  i.e., 
individual  E had  disappeared  from  the  observer  at  time  6.5  (so  that 

the  fact  of  E's  deatli  at  time  7 is  unknown)  . 

^ P2(t)  - Modified  data 


0-1 

1.0 

1-3 

0.8 

3-6.5 

0.533 

Since  tlie  time  of  the  death  for  individual  E is  now  unknown, 
one  can  only  estimate  that: 

0 £ P2^^^  1 0.533  for  t > 6.5 
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If  the  analyst  is  willing  to  assume  a functional  form  for  the 
survival  function,  then  he  may  calculate  the  manner  in  which  the 
estimator  P2(t)  decreases  to  zero.  However,  the  data  alone  are  insuf- 
ficient when  a strictly  non-par ame trie  estimator  is  used. 

The  product-limit  estimator  is  a useful  and  intuitively  appeal- 
ing method  of  dealing  with  incomplete  observations.  It  has  been  wider 
used  and  studied.  However,  the  product-limit  has  one  disturbing 
characteristic ; 

Most  of  the  biological,  physical  or  other  causes  of  deaths  pro- 
duce a survival  probability  that  continuously  decreases  in  time. 

It  is,  therefore,  one  may  be  a little  uncomfortable  estimating 
the  survival  probability  with  a step  function.  One  is  tempted 
to  smooth  the  estimator  to  make  it  a monotonic  decreasing  func- 
tion of  t. 

3.  The  Third  Estimator,  "p^Ct)" 

p^it)  is  a modification  of  p^Ct).  Like  P2(t),  it  is  a step 
function  with  discrete  drops  at  those  times  corresponding  to  the  observed 
deaths  in  the  sample  population.  It  may  also  be  expressed  as  a product 
of  conditional  probabilities: 

i 

P3(ti)  =TrP3(\lVl)  (E-4) 

3=1 


where  the  t are  the  times  of  observed  deaths  and  t is  zero.  The  con- 

K O 

ditional  probabilities  on  the  right-hand  side  of  Equation  (E-4)  differ 
somewhat  frcrni  those  in  Equation  (E-2) : 


N^-1 


(E-5) 
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Equation  (E-5)  differs  from  Equation  (E-3)  in  the  interpreta- 
tion of  the  numbers  of  individuals  at  risk.  Here,  the  value  of  N,  is 

k 

taken  to  be  the  average  number  of  individuals  observed  surviving  in 
the  interval  between  the  (k-l)£t  observed  death  and  the  kth  observed 
death.  The  number  of  observed  survivors  decrease  at  intermediate  times 
if  events  are  censored,  and  hence  the  are  not  necessary  integers. 

The  value  of  is  regarded  as  the  effective  sample  size  for 
the  interval  from  tj^  ^ to  t^^.  In  the  sample  data  base  shown  above, 
individual  B is  known  to  have  survived  from  time  1 to  time  2,  or  half 
of  the  interval  between  the  first  death  at  t=l  and  the  second  death  at 
t=3.  Therefore,  the  estimator  p^  treats  individual  B as  half  a parti- 
cipant in  the  interval  between  the  death  of  individuals  A and  C. 


The  effective  sample  size  for  this  interval  is  then  3.5 
(2-1) 

(n  = 3 + ---  = 3.5)  (full  contributions  from  individuals  C,  D and  E, 

plus  a half  contribution  from  B) . For  our  hypothetical  data  base,  the 


following  values  are  calculated  for  p. 


P3(t) 


5/5  = 1.0 


1-3  4.5  X 1.0  = 0.8 

3-7  (2. 5/3. 5)  X 0.8  = 0.571 

(7)  (1.75/2.75)  X 0.571  = 0.364 

The  value  of  P3(t)  can  never  be  less  than  the  corresponding  value  of 
P2(t).  In  the  special  case  with  no  censoring  events  the  estimators 
Pj^(t),  P2(t)  and  P3(t)  are  identical. 

One  might  perturb  the  data  by  shifting  the  time  of  B's  trunca- 
tion event  down  to  1+e  or  up  to  3-e,  e arbitrarily  small.  The  depend- 
ence of  the  estimator  p^  upon  the  exact  time  of  the  censoring  events 
may  now  be  demonstrated. 
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For  purposes  of  illustration,  the  tine  of  the  censoring  event 


for  individual  B(t2)  is  decreased  from  2 to  1.1,  then  increased  to  2.9. 


t 

P3(t),  = 2 

P3(t),  t^  = 1.1 

U) 

ft 

! ft 
1^ 

0-1 

1.0 

1.0 

1.0 

1-3 

0.80 

0.80 

0.80 

3-7 

0.571 

0.538 

0.597 

(7) 

0.364 

0.342 

0.380 

This  example  demonstrates  an  intuitively  appealing  characteris- 
tic of  the  estimator,  p^.  As  the  total  observed  survival  time  increases 
for  the  individuals  in  our  sample  (with  deaths  held  constant) , the  value 
of  the  estimating  function  increases  over  at  least  a portion  of  its 
range. 

We  may  safely  assume  that  the  true  survival  function  eventually 
tends  to  zero  with  time,  since  no  physical  or  biological  system  lives 
forever.  However,  there  are  no  observations  on  the  survival  of  indi- 
viduals beyond  time  7.  The  data  only  indicate  that  our  step-function 
estimator  drops  to  a value  of  .364  at  t=7,  but  the  nonparametric  esti- 
mator gives  no  information  about  the  survival  function's  subsequent 
decline  from  .364  to  zero.  However,  the  data  alone  are  insufficient 
when  a strictly  nonparametric  estimator  is  used. 

B.  POINT  ESTIMATOR 

As  mentioned  above,  the  estimators  p^^,  P2  and  p^  are  somewhat  unde- 
sirable because  they  give  step-function  estimates  for  a continuous 
survival  function.  The  next  three  estimators  p^,  p_  and  p^  are  modifi- 
cation  of  the  first  three.  Again  they  provide  estimates  of  the  survival 
function  only  at  those  points  in  time  that  corresponde  to  observed  deaths. 
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These  estimators  are  specified  by  Equations  (E-2)  and  (E-4) , except  for 
a substitution  of  the  term  (N+1)  in  place  of  (N) » 


Since  the  point  estimators  have  rigorous  definitions  at  only  dis- 
crete points  in  time,  it  is  necessary  to  offer  an  interpolation  rule. 

That  is,  we  need  a method  of  "connecting  the  dots."  The  method  proposed 
here  is  to  assume  that  the  survival  function  declines  in  a piece-wise 
exponential  decay  between  the  discrete  points  in  time.  This  procedure 
is  equivalent  to  assuming  that  the  hazard  function  is  essentially  con- 
stant between  a consecutive  pair  of  the  discrete  times,  but  that  the 
hazard  varies  from  one  time  period  to  the  next.  Such  an  assumption  is 
intuitively  acceptable  unless  one  suspects  violent  fluctuations  in  the 
hazard  function. 

1.  The  Estimator,  "p^(t)" 

p^(t)  is  analogous  to  p^(t)  in  that  only  those  individuals  ob- 
served to  die  are  included  in  the  sample.  These  two  estimators  are  naive 
because  they  suppress  all  data  from  the  survival  times  of  individuals 
terminated  from  observation  by  censoring. 

These  estimates,  i.e.,  p^(t)  and  p^(t),  tend  to  ignore  informa- 
tion from  the  more  long-lived  individuals  in  the  sample,  and  they  may 
be  expected  to  give  biased  estimates  of  the  survival  function. 

The  point  estimator  p^(t)  gives  the  following  values  with  sample 
data  base  presented  earlier  in  this  section. 


t 

ft 

Interpolation 

t 

0 

— 

1.0 

f^n(3/4) 

e 

1 

3/4  = 0.75 

■ 

^ •i.n{2/2) 

3 

(j)  X 0.75  = 0.5 

1 

. ^ •Zn(l/2) 

7 

(-)  X 0.5  - 0.25 

l__£ 

Hi 

(i)e 
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I 


A 


2.  The  Estimator,  "p^(t)" 

D 

The  estimator  Pg(t)  similarly  corresponds  to  the  product-limit 
estimator  p^Ct).  These  two  estimators  use  information  from  the  indi- 
viduals on  whom  there  are  censored  observations,  p^ , like  p^,  does 
not  exploit  information  about  that  portion  of  the  censored  observation 
after  the  death  event  (of  some  other  individual)  preceding  the  censor- 
ing event. 
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For  our  hypothetical  data  base  the  following  values  are 
calculated  for  p^(t): 


V'fhenever  censored  observations  are  present,  the  estimator  p^(t)  never 
exceeds  p^ (t) . 

For  p^(t),  the  value  of  is  taken  to  be  the  number  of  surviv- 

ing individuals  in  the  saraple  just  before  the  observation  of  the  ith 
death.  This  value  is  smaller  than  the  number  of  surviving  individuals 
just  after  the  (i-l)s_t  death  if  any  truncation  events  occur  in  the 
interval.  In  fact,  is  the  smallest  number  of  surviving  individuals 
observed  at  any  time  during  the  interval  (t^  t^) . Thus  p^  might  be 

expected  to  introduce  a bias  by  using  values  of  that  are,  on  the 
average,  too  small.  However,  this  bias  would  be  much  less  severe  than 
the  bias  anticipated  for  the  estimator  p^(t) . 

The  estimators  and  p^  are  insensitive  to  the  precise  times 
of  the  censoring  events.  A change  in  the  time  of  the  censoring  event 
for  individual  B to  l+£  to  3-e,  e arbitrarily  small,  does  not  alter  the 
estimates  from  p^  and  p^  given  in  the  preceding  paragraph. 
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3.  The  Estimator,  "p^(t)" 

The  estimator  (t)  corresponds  to  p, (t)  by  accounting  for  all 

D J 

of  the  survival  time  for  the  truncated  observations.  For  our  hypo- 
thetical data  base,  the  following  values  are  calculated  for  p- (t) ; 

6 


The  estimator  p (t)  is  based  on  the  average  number  of  surviving 

individuals  noted  in  the  various  time  intervals.  These  estimators  give 

part  credit  for  individuals  whose  lifetime  is  censored  in  mid-interval. 

The  value  of  tl.  for  p, (t)  is  an  unweighted  time  average.  If  the  obser- 
1 b 

vation  of  an  individual  is  truncated  after  23%  of  the  interval  has 
elapsed,  then  that  individual  contributes  a value  of  0.23  to  N^.  Indivi- 
duals who  are  observed  to  survive  the  entire  interval,  and  the  individual 

whose  death  terminates  the  interval  each  contribute  a value  of  1.0  to  N.. 

1 

This  interpretation  of  the  effective  sample  size  is  approximate  if  the 
hazeurd  is  approximately  constant  over  the  interval.  If  the  hazard  function 
changes  markedly  within  a time  interval  containing  censored  events,  then 
this  interpretation  of  the  effective  sample  size  is  biased.  Therefore, 
the  procedure  of  determining  the  value  of  tl.  for  the  estimator  p,  (t)  is 
based  on  the  implicit  assumption  that  the  survival  function  is  locally 
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l' 


[ 


I 


exponentialo  If  the  hazard  function  may  be  assumed  to  vary  slowly  over 
each  of  the  time  intervals  (t.  , t.)  then  p would  appear  to  be  biased 

1”X  1 D 

on  an  acceptaible  approximation. 

The  estimator  p , like  p,,  depends  on  the  precise  times  of  all 

O J 

deaths  and  censoring  events. 


P6(t),  = 2 


Pg  (t) , t^  = 1.1 


Pg(t),  t^  = 2.9 


0 


1.0 


1.0 


1.0 


1 5/6  = 0.833 
3 (-j^)  X 0.833  = 0.648. 
7 0.648  = 0.412. 


5/6  = 0.833 

(|^)  X 0.833  = 0.628. 
(■i^)  X 0.628  = 0.399. 


5/6  = 0.833 

3 95 

(j^)  X 0.833  = 0.665 
X 0.665  = 0.423 


This  illustrates  that  an  increase  (or  decrease)  in  the  total  observed 
survival  time  causes  an  increase  (or  decrease)  in  the  estimate  p^  over 
at  least  some  of  its  time  range. 

If  the  last  event  is  a censored,  and  not  an  observed,  death, 
these  estimators  also  require  definition  for  the  time  period  starting 
with  the  time  of  the  last  death  and  ending  with  the  time  of  the  final 
censoring  event. 

The  method  proposed  here  for  P^(t)  and  Pg(t)  is  to  continue  the 
exponential  function  used  in  the  interval  terminated  by  the  time  of  the 
last  death.  This  procedure  can  be  illustrated  with  the  modified  data 
base  used  above  in  the  discussion  of  P2  and  p^. 


C.  THE  BAYESIAN  ESTIMATORS 

Consideration  is  next  given  to  quasi-Bayesian  estimators  based  on 

a uniform  prior  distribution  on  the  unit  interval.  Let  X ,...,X  be  the 

1 N 

true  survival  times  of  N individuals  which  are  censored  on  the  right  by 
N follow-up  times  Y^,...,Y,j.  It  is  assumed  that  the  X^  are  independent. 
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identically  distributed  random  variables  with  common  distribution  p(t) 
and  we  wish  to  estimate  the  survival  function 


p(t)  = Pr(x  > t) 


However,  we  only  have  available  the  data. 


Z.  = min  (x. , Y. } 
1 11 


1 if  X.  < Y. 

1 — 1 


0 if  X.  > Y. , i=l, . . . ,n 
1 1 


If  6.  =0,  then  Z.  is  called  "a  loss",  and  if 
1 1 

6.  =1,  then  Z.  is  called  "a  deatli". 

1 1 

Then  p [6.  =1]  = p [X.  > t]  = p(t),  i=l,...,M- 
r 1 r 1 

The  maximum  likelihood  estimator  for  p(t)  is 


p(t)  = — where  s = E 6. 

” i=l  ^ 


is  the  number  of  successful  tests,  s has  the  binomial  distribution. 


P{s|p)  — p (l“p,*  , s=0 , 1 , . . . , 0 < p < 1 

fp(p)  = 1,  0 < P < 1 


The  joint  density  of  s and  p is 


^s,p^®'P^  = (^)  p®{l-p)“"®,  0 < p < 1,  s=0,l,...  N. 


The  marginal  for  s is 


s..  . n-s , . s!(N-s)!  1 

Pg(s)  = f p (1-p)  dp  = (g)-  (nVDT  ■ nTI 
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for  s=0,l,...N.  Thus,  averaging  over  the  values  of  p,  all  of  which  are 
assumed  to  be  equally  likely,,  ^he  values  of  s are  equally  likely  to  occur. 
The  posterior  for  p then  is 

I ^ (tI+2)  s U-s 

fp|3(pls)  =r(T+i)f('N-s+i)  P , 0 < p < 1 , 

a beta  density  witli  parameters  s+1  and  N-s+1.  The  mean  of  the  posterior 
is  (s+l)|  (^J+2)  and  the  modal  (maximum  value)  of  the  posterior  is  s/li ; thus 
the  Bayes  estimate  of  p (given  s survivers  occur  in  the  sample  of  N ) is 


* s+1 

P “ H+2 


(C-1) 


Then,  equation  (C-1)  yields  a step  function  and  also  has  shown  that  the 
uniform  prior  has  the  effect  of  adding  two  individuals  to  the  popula- 
tion at  risk  with  one  dying  at  time  zero  and  the  other  essentially 
immortal. 

The  Bayesian  estimators  based  on  a uniform  prior  distribution  on 
the  unit  interval  are  denoted  E^j^(t)  , and  |^^{t)  , that  correspond, 

respectively,  to  the  estimators  p^ (t) , P2(t)  and  P2(t).  The  sample 
data  base  thus  gives  the  following  estimates  of  the  survival  function: 

t Pll^*^^  Pl2<^>  Pl3<^> 


0-1 

4/5  = 0.8 

6/7  = 0.857 

6/7  = 0.857 

1-3 

3/5  = 0.6 

X 

0.857  = 0.714 

X 

0.857  = 0.714 

3-7 

2/5  = 0.4 

X 

0.714  = 0.536 

X 

0.714  = 0.556 

(7) 

1/5  = 0.2 

X 

0.536  = 0.268 

1.75 

^2.75 

)x 

0.556  = 0.354 
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At  the  time  of  the  final  event  (whether  a death  or  a truncation) 


these  step-function  estimators  drop  to  some  positive  value.  Again, 
we  have  no  data  to  indicate  how  the  survival  function  proceeds  to  zero 
at  subsequent  times. 


I 


D.  THE  JACKKNIFE  ESTIMATOR 


I 


We  will  assume  that  we  observed,  or  have  generated  in  a simulation, 
a survival  probagility  p(tj),  j=l,...,n,  from  various  sample  sizes. 
Furthermore  we  have  some  parameter  or  characteristic  plt^)  of  the 
sample  size  which  we  wish  to  estimate  with  an  estimator  p(tj).  The 
jackknife  estimator  p{t,n)  described  below  is  an  approximately  unbiased 
estimator  of  pCt^).  A modification  of  it  has  other  useful  properties. 

p_^(t,n-l)  is  the  estimator  from  the  sample  of  n of  the  X^'s  with  the 
ith  value  deleted  from  the  sample. 

p^(t,n)  = n p(t,n)  - (n-1)  p j^(t,rv-l)  i=l,...,n 

p(t,n)  = — E p.  (t,n)  = n p(t,n)  - — E p , (t,n-l) 

n . , 1 n . , -i 

1=1  1=1 

the  p^(t,n),  called  the  PSEUDO-values . 

The  PSEUDO-values  can  be  used  to  obtain  variance  estimates  of  p(t,n) 
and  to  set  approximate  confidence  limits,  using  Student's  t. 

The  idea  is  that  the  PSEUDO-values  will  be  approximately  indepen- 
dently and  normally  distributed.  The  jackknife  estimator  p(t,n)  is 

2 

sample  average  so  we  form  an  estimate  S,s  of  its  variance  given  by 

p(t,n) 

the  following  relationship  (Miller,  1974) : 


^p(t,n) 


Ep^^(t,n)  - i(Ep^ (t,n) ) ^ 
n-1 


n 
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L 


This  procedure  is  particularly  useful  if  the  number  of  data  points 

is  small,  but  it  must  be  used  with  care.  Note,  that  the  estimator  p(t,n) 

1 

is  designed  to  eliminate  a — bias  term  in  the  estimator  p(t,n).  Of 
course  the  computational  aspects  of  the  complete  jackknife  can  be  quite 
onerous,  especially  if  p(n)  were,  say,  a complicated  maximum  likelihood 
estimator.  Miller,  reference  (4)  has  shown  that  the  product  limit 
estimator  is  its  own  jackknife. 

Logistic  Transformation 

Although  one  can  legitimately  jackknife  the  Kaplan-Meier  estimate 
directly,  there  is  some  reason  to  believe  that  a preliminary  transforma- 
tion will  give  improved  results.  Consequently,  consider  the  transforma- 
tion 


a = In  { 


. pit) . . 

l-p{t) 


) 


and  notice  that  where  the  range  of  p(t)  is  from  zero  to  unity,  the  above 
transformation  makes  the  range  of  Z run  from  -<»  to  The  procedure 
utilized  will  be  as  follows. 

(A)  Compute  the  overall  estimate  at  a time  point  t,  using  all  N data 
points,  and  using  a “continuity"  correction  that  has  the  effect 
of  removing  the  effect  of  a zero  in  the  logarithm  (see  D.R.  Cox, 
Analysis  of  Binary  Data,  Methuen  Monograph) ; 


Z = Zn( 

N 


-) 


(B)  Compute  the  ^.-values  by  leaving  out  each  data  point  in  turn  when 
computing  p(t):  for  i=l,2,...,N. 


Z = iln('  2(H;12 ) 
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2(N-1) 


(C)  Form  the  pseudo-values 


2.  = NJl„  - (N-1)  £ , 

1 N N-1,-1 


(D)  Compute  z,  S 

z 

(E)  Put  approximate  confidence  (l-a)*100%  limits  on  E[Jl]  as  follows 


where 


L < E [£]  < H 

H(L)  = z +(-)  t^_^  (N-1)  / ^ 


(F)  Transform  bash  to  obtain 


L H 

e , e 

T-fT- ' 

1+e  1+e 


The  true  value,  p(t),  should  be  enclosed  between  these  levels  for 
roughly  (l-a)*100%  of  all  samples.  The  coverage  properties  of  this  pro- 
cedure will  now  be  checked  by  simulation;  successive  samples  of  size  N 

will  be  selected,  the  jackknife  limits  H and  L will  be  computed  for  each, 

L H 

6 @ 

and  a check  will  be  made  as  to  whether  — < p(t)  < — or  not.  Tables 

- L ” ” _ H 

1+e  1+e 

illustrating  performance  are  given  subsequently. 


Let 


il.,  , . is  the  logistic  transformation  estimator  from  the 

N-1,-1 

of  the  X^'s  with  the  i;^  value  deleted  from  the  sample. 


1 

N-1 ,-L 


^N-l,-i  2 (N-1)  ) 


^“^N-l,-i^^^'^2(tJ-l) 


N-l,-i 


i 

t 

1 

2 

3 

4 

5 

3.04 

0.98 

0.98 

0.98 

0.98 

^2 

3.04 

0.98 

0.98 

0.98 

0.98 

0.63 

0 

0.98 

-0.46 

-0.46 

^4 

0.63 

0 

0.98 

-0.46 

-0.46 

-3.04 

-3.04 

-3.04 

-1.89 

z.  = Nt  - (N-1)  I . 

X N N-1,-1 


= N2.n  (■ 


PN<t)  + 2n 


2N 


-p)  - (N-1)  IL^( 


Pn-I,-!^^^  ^ 2(M-1)  ^ 

^“^N-l,-!^’^^'^  2 (N-1) 


z^(U)  are  called  PSEUDO-values  of  logistic  transformation, 


sample  n 


the 


following  values  are  calculated: 


Average  of  the  pseudo-values 


- 1 

z = — E z. 
i=l  " 


Invert  to  find  jackknife  estimator  of  logistic  transformation 


p (t) 


= Jin 

l-p(t).- 

1 , z — 

_ ( 1-*-zTtJ  e - 2U 

1 + e ^ 


called  the  jackknife  estimator 
of  logistic  transformation 


Variance  of  the  z 


S = Var( z ) 

z 


1 

n-1 


n 

E 

i=l 


Zi  - z 


The  following  values  are  calculated: 


B 

z 

Var  J 

B 

0,5484 

0.646 

13.6 

"2 

0.54Q4 

0.646 

13.6 

"3 

0.0568 

0.516 

6.727 

"4 

0.0568 

0.516 

6.727 

-3.882 

0 

3.361 

The  jackknife  estimator  for  estimating  variability  and  giving  confidence 
interval . 

Tukey,  reference  C3)  has  suggested  that  in  the  jackknife  procedure 
wo  consider  tJie  pseudo  values  z^(n)  as  approximately  independent  and 
identically  dictricuted  and  consequently,  since  z is  an  average  of 
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r 


the  Z ^ (N) , proceed  as  if 


z - Z 


i=l 


has  t-distribution  with  N-1  d.t’ . 

If  the  are  approximately  normal  variates  (Miller  has  shown) 
confidence  bands  for  the  unknown  p(t)  are  given,  as  for  the  mean  of 
any  normal  variate  when  estimated  from  scunple  size  n. 


" i 'l-a/2 


i.e. 


_ s p (t)  + - — s 


1-P(t)+  - 


L(n)  - z =■  t,  „ 

^ l-“/2 


L(n)  = z + — t 


l-a/2 


(1+  - A.  A.,eL(N)  _ 

' < p,t)  < 211 


1 + ei®’  1 + 


The  following  values  are  calculated; 


4 

Va/2  = 

Lower  Int. 

Upper  Int. 

"l 

0 

1.0 

0 

1.0 

0 

1.0 

% 

0 

1.0 

S 

0 

0.14 
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The  basis  for  this  leap  of  the  imagination  seems  to  be  that  if 
^ = X = then  the  procedure  for  obtaining  confidence  intervals 
using  equation  (D-1)  and  pseudo-values  is  the  same  as  the  procedure 
using  jackknife.  Then  if  % = ^ and 


z 


z . 
1 


we  have 


z , 
1 


= N£  - 
N 


(n-1)  £ 

N-l,-i 


= N X.  - '(N- 1 ) 


A 
{ E 
j=l 


X.  } - X, 


N-1 


N 

= E 
j=l 


X . 
D 


N 

[ E X.]  + X,  = X, 

. 1 1 1 1 
3=1 


Thus  the  pseudo  value 

1 

z.  = X.  and  z = - E X.  = X 
1 n . , 1 n 

1=1 

The  pseudo  values  are  independent  if  z = X^  and  they  are  normal  if 

X.  is  normal. 

1 

E.  PARAMETRIC  ESTIMATOR,  "p^{t)" 

This  paper  considers  one  additional  estimator,  denoted  p^  (t) . It 
is  a parametric  estimator.  Therefore,  it  is  not  really  a competitor  to 
the  thirteen  non-parametric  estimators  considered  here.  In  general,  a 
parametric  estimator  would  not  be  used  if  the  functional  form  were  re- 
garded as  unknown.  Similarly,  a non-parametric  estimator  would  not 
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normally  be  used  if  the  survival  function  were  strongly  suspected  to 
have  a specified  form. 

p^(t)  is  the  well  known  maximum  likelihood  estimator  for  the  expo- 
nential distribution; 


where 


p.^(t)  = e 


-t/T 


I t. 

1 

number  of  observed  death 


In  our  sample  data  base,  the  total  observed  survival  time  is  19,  and 
three  deaths  are  observed.  Thus , 


and 


Zt^=l+2+3+6+7=19 


-3t/19 


Calculations  for  selected  times  of  interest  yield  the  following  esti- 
mates : 


p^(0)  = 1.0 
p^(l)  = 0.854 
p.^(3)  = 0.623 
p^(7)  = 0.331 

The  thirteen  non-parametric  estimators  are  compared  for  a variety  of 
generating  distributions  for  both  the  death  mechcinism  and  censoring 
mechanism. 
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IV.  INSTRUCTIONS  FOR  USING  PROGRAM 


f INPUT 

I 

Each  input  card  bears  nine  variables.  The  distribution  of  time  of 
death  is  entered  in  the  first  set  of  (five)  columns,  the  censoring  dis- 
tribution is  entered  in  the  second  set  of  (ten)  columns,  a parameter  of 
the  censoring  distribution  is  entered  in  the  third  set  of  (ten)  columns, 
the  number  of  replication  is  entered  in  the  fourth  set  of  (five)  columns, 
the  number  of  the  event  is  entered  in  the  fifth  set  of  (five)  columns. 
For  the  purpose  of  all  print  output  used  code  "0"  and  "1"  in  the  sixth 
set  of  (five)  columns,  the  seed  number  is  entered  in  the  seventh  set  of 
(five)  columns,  after  the  card  giving  the  time  of  the  last  event  of  a 
data  set,  a card  with  "0"  or  "1"  in  the  column  50  is  inserted,  i.e.,  the 
"0"  indicating  more  data  sets  to  follow  and  "1"  indicating  the  last  data 
sets  and  t value  is  entered  in  the  ninth  set  of  (eight)  columns. 

The  distribution  of  timeof  death  and  of  censoring  time  used  code  as 
follows; 

Code  Type  of  Distribution 

1 Uniform 

2 Exponential 

3 Delta  function 

OUTPUT 

The  output  lists: 

1)  the  time  of  each  observed  failure 

2)  estimated  survival  probability  at  that  time 

3)  the  variance  of  that  estimator 

4)  result  of  goodness  fit 


L 
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a)  mean  error 


I 


b)  mean  absolute  error  (ABS) 

c)  root- mean- square  error  (RMS) 

5)  total  number  of  observed  death 

6)  confidence  interval  at  particular  time 
Definition  of  Fortran  Variables 

NDIE  : the  distribution  of  time  of  death 

NTRUNC  : the  distribution  of  censoring  time 

XTRUNC  : the  parameter  of  the  distribution  of  censoring  time 

NREPL  : niamber  of  replication 

NEVENT  : number  of  event 

NWRITE  : write  all  output  or  partial  output  of  simulation 

NEND  : indicate  more  data  sets  or  last  data  set 

TN  : t statistic  value 

Pj^  : the  estimator,  Pj^(t) 

p^  : the  estimator,  p^ (t) 

p^  : the  estimator,  p^ (t) 

p^  : the  estimator,  p^(t) 

p^  : the  estimator,  p^(t) 

p.  : the  estimator,  p, (t) 

6 6 

p^  : parcimetric  estimator  , p^{t) 

p : jackJcnife  estimator  of  logistic  transformation  of  p,(t) 

Pg  : jackknife  estimator  of  logistic  transformation  of  p^Ct) 

p.  : jack)cnife  estimator  of  logistic  transformation  of  p,.  (t) 

p^j^  : Bayesicin  estimator  of  p^(t) 

Pj^2  = Bayesian  estimator  of  P2  (t) 

p^2  : Bayesian  estimator  of 


39 


: jackknife  estimator  of  logistic  transformation  ofp^Ct) 

SL'T,J)  ; PSEUDO- value 

SBAr  ; average  of  pseudo-value 

Var  : variance  of  estimator,  p(t) 

Var  : variance  of  jackknife  estimator 

J 

u(I,J)  : mean  of  goodness  fit 

w(I,J)  : edasolute  mean  of  goodness  fit 

s(I,J)  : root  mean  square  error 

: upper  confidence  interval  of 

: lower  confidence  interval  of 

C : upper  confidence  interval  of  p_(t) 

3 o 

: lower  confidence  interval  of  Pg(t) 

: upper  confidence  interval  of  Pg(t) 

Cg  : lower  confidence  interval  of  Pg(t) 

: upper  confidence  interval  of  p^^^  (t) 

Cq  : lower  confidence  interval  of  p^^Ct) 

o 10 
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To  compare  RMS  with  product  limit  (p  (t))  and  jackknife  estimator  of 

logistic  transformation  (p, .(t)) 
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Computer  output  of  the  fourteen  estimators 
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V.  RESULTS  OF  THE  SIMULATION 


This  paragraph  presents  graphical  comparisons  of  the  13  estimators 
based  on  simulated  data.  Each  comparison  is  based  on  1000  replications 
of  a simulated  data  base.  The  bias  and  RMS  error  (square-root  of  mean- 
squared  error)  of  each  estimator  depends  on  the  parameters  that  control 
the  simulated  data  base.  No  single  estimator  dominates  all  others 
under  all  conditions. 

The  bias  and  RilS  errors  of  the  estimators  depend  on  several  factors; 

(A)  The  sample  size  (NEVENT)  of  individuals  under  observation  at  time 
zero  affects  the  accuracy  of  the  estimators.  In  general,  a larger  sample 
size  leads  to  a better  estimate  than  a smaller  sample.  Values  of  NEVENT 
selected  for  simulation  are  5,  10,  25,  and  50  (plus  one  simulation  with 
NEV’ENT  = 100)  . 

(B)  The  distribution  of  times  at  which  the  observations  are  censored 
(unless  the  individual  dies  earlier)  affects  the  performance  of  the 
various  estimators.  This  distribution  is  particularly  important  in  con- 
junction with  the  distribution  of  lifetimes  (do  most  individuals  die 
before  censoring  is  li)cely?,  are  deaths  and  censoring  events  about 
equally  likely  at  all  times?,  are  most  observations  censored  before  death?). 
Three  types  of  distributions  are  assumed  to  underlie  the  censoring  mech- 
anism: 

(1)  Some  of  the  samples  are  generated  on  the  assumption  that  no  censor- 
ing occurs. 

(2)  Some  samples  are  generated  from  a uniform  distribution  of  times  of 
censoring. 
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(3)  Other  data  bases  are  generated  from  an  exponential  distribution  of 


censoring  times. 

(C)  The  distribution  of  lifetimes  (ignoring  the  possibility  of  censor- 
ing) also  affects  the  performance  of  the  various  estimators.  Two  types 
of  distributions  are  assumed  to  underlie  the  death  mechanism; 

(1)  Some  of  the  samples  are  generated  from  a uniform  distribution  of 
lifetimes. 

(2)  Other  data  bases  are  generated  from  an  exponential  distribution  of 
lifetimes . 

If  a uniform  distribution  of  lifetimes  is  selected,  its  range  is 
always  over  the  interval  from  time  0 to  time  1.  If  an  exponential  dis- 
tribution is  selected,  it  always  has  a mean  lifetime  of  1.  The  distri- 
butions of  truncation  times  (uniform  or  exponential)  have  parameters 
.25,  .5,  .667,  .75,  1,  1.333,  1.5,  2 and  4.  A wide  variety  of  seutiples 
may  be  simulated  by  mxing  various  pairs  of  distributions  (for  censoring 
times  and  deaths) . Since  the  time  units  are  arbitrary,  the  restriction 
on  mean  lifetimes  is  irrelevant. 

The  true  value  of  the  survival  function  is,  p(t),  and  the  form  of 
this  function  affects  the  relative  performance  of  the  13  nonparametric 
estimators.  For  example,  the  Bayesian  estimator  tends  to  be 

better  as  measured  by  square-root  of  mean-squared  error  than  its  counter- 
part (the  product-limit  estimator,  P2(t))  for  the  time  frame  in  which 

.3  < p(t)  < .9 

However,  the  product  limit  estimator  tends  to  be  better  for  those  times 
when  p(t)  is  close  to  zero  or  unity. 

The  point  estimators,  Pc(t)  and  p, (t)  tend  to  be  better  than  the 

D O 

product-limit  estimator  (P2(t))  for  all  time  periods.  The  jac)c)cnife 
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estimators  of  logistic  transformation  (Pg(t),  Pg(t),  Pj^Q^t))  of  point 

estimators  tends  as  same  as  its  counterpart  point  estimators  (p^(t), 

p (t) , p (t) ) for  all  time  periods.  And  also  the  estimator  formed  by 
b 6 

jackknifing  the  logistic  transformation  (Pj^^(t))  of  the  product  limit 
estimator  tends  to  be  better  than  its  counterpart  product  limit  (P2(t)) 
for  the  time  frame  in  which 

.1  < p(t)  < .7 

However,  the  product  limit  estimator  tends  to  be  better  for  those  times 

when  p(t)  is  close  to  unity.  Point  estimators,  P-Ct)  and  p^  (t)  tend  to 

b 6 

be  same  for  the  time  freime  in  which 

0.1  < p^(t)<  0.9 

However,  the  P2{t)  tends  to  be  better  for  those  times  when  p^(t)  is 
close  to  unity.  The  jackknife  procedure  may  be  validated,  in  an  empiri- 
cal sense,  by  sampling  experiments  or  computer  simulation  in  the  follow- 
ign  manner.  First,  times  of  censoring  and  death  are  obtained  by  drawing 
random  numbers  from  postulated  distributions.  Second,  the  jackknifed 
estimator  of  the  logistic-transformed  product-limit  estimation  is  found, 
and  confidence  limits  are  computed  by  the  method  of  Tukey,  reference  (3) . 
Since  the  true  value  of  survival  function,  p(t),  is  known,  so  is  the 
theoretical  value  of  A.  The  jackknife  confidence  intervals  can  be 
checked  for  coverage;  if  £ A ^ then  the  particular  interval  covers, 
while  otherwise  (if  A < or  < A)  it  does  not  cover.  Finally,  the 
above  procedure  can  be  repeated  many  times  (say  1000)  and  the  fraction  of 
repetitions  which  contains  the  true  value  of  A is  recorded.  This  fraction 
of  the  coverage  should  desirably  be  close  to  (1-a) , the  nominal  confidence 
level. 
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The  jackknife  confidence  limit  procedure  can  be  said  to  be  robust 
of  validity,  ref  (7),  if  the  actual  coverage  is  close  to  the  nominal 
coverage,  1-a,  for  a various  distributions.  Such  seems  to  be  true  for 
large  n (n  ^ 50) . However,  the  jackknife  confidence  limits  do  not  cover 
accurately  when  the  true  value  of  p(t)  is  close  to  unity. 

The  following  tables  illustrate  confidence  limits  of  jackknife 
method  of  product  limit  (p^Ct)) . Many  computer  generated  graphics  are 
presented  on  the  following  pages  to  complete  this  section. 
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51 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FIMISHED  TO  DOC  


• ( t i 


. / » • ; • s . j 


. 'Jr  ‘ 


I . 7 / ' 


♦ l.'77r-ci 


i. 


l.tl  Cf--01 


1. 


I . <ir7r-.j! 


l.2-< 


ur^t-L-oi 


I . )o  - )l 


u. ' < 


• i 1-  <-  •<  f V t* 
0.  7S 


l.-.^ 


\.Cf  cr-ci 


Fig.  2;  Comparison  of  root  mean  square  derived  from  sample  size  5 
(Step  function  estimators  vs.  jackknife  estimators  of 
logistic  transformation) 


52 


t 


‘vO'  1*  I ".fM 


r; i . J . v^.  r;  t 


..'.:.:.i  \5vrMT^ 


t ) . * ,*  1 .'  1 ■ ) . V'* , ‘‘  1 3t  * J 


^..,f  • f**  t V « v'.  ♦<  ^***  < V Cl.*  V.  V 


z.t  C*‘'I  -vU 


l.J  -Ol 


l.37K-n  ♦ 


I ♦ 

• ^ 4 • < • ' 

O.T 


l.OJ 

i > i.  V •'  « <“ 


I 4 l t t i '.■■  f 

•»,  n.so 


l.OO 


.3.f » ‘^r-c! 


7,6CK-Ol 


1.  Ff  ^‘■-ci 


1.  -vl 


1.  inF-oi 


Fig.  3;  Comparison  of  root  mean  square  derived  from  sample  sizes. 
(Step  function  estimators  vs.  Bayesian  estimators) 


53 


THIS  PAGE  IS  best  QUAIITY  FRACIICABLZ 
RaJM  COPY  FURNISHED  TO  DDC  ^ 


n ■ i /♦ . u .»«•  Nf  v<  - ID 


'*  •V‘>« 


I -0! 


->  • t i I • I '.) 

♦ 


I . - D) 


1.  I33t- Jl 


1.  i-ar-ci 


I. 


♦ i.26{,r-ji 


1. 


^ . j oof  - Jr 


J • c- 


1. 30 


f .6<fC-02 


Fig.  4:  of  root  mean  square  derived  from  sample  size  10. 

(Step  function  estimators  vs.  point  estimators) 


PT  IS  BEST  HUAETTE 


I M 1 i»  Uf  . - ’ XI.  l..xl  - 


»‘.'  I • » . V >.  .’  ..  i » J . V ..  r‘  . ( O , V . . 01  ) ( X > 


.f  \ = 1 J 


I .*>  701  -J\  ♦ 


I . 0-  - 


I - <»  1 J I - M t 


1 • 1 jii  • - I 


1 • J4cr  - Ji 


1 1 ^ K f<  f*** 


• 1.6701-01 


K ♦ 

o.'..) 


« « i;  f F 


6.  f'i 


1 .(Ht 


Fig.  5;  Comparison  of  root  mean  square  derived  from  sample  size  10, 
(Step  function  estimators  vs.  jackknife  estimators  of 
logistic  transformation) 


♦ l.-ilOF-Ul 


« -I'l 


e.^-cci  -0^ 


55 


r 


S CO.?  mcnoAwi. 

WUM  COf  X i URNlSHiD  10  DOC 


N .1’  I * I » I r‘  i M 


M V i f J ^ 10 


^4  I . I . ^ >.  * » ' » • ' 


>bi  • 'I 


I . ? f -ul 


I . ; t *1  - J ! 


l . l J3I  -01 


\ . ^0 

i < »•  t ♦ < A 


l.c( * » -Jl 


I . * < ^ ‘ : I 


1 . 1 ' ■»  f - 0 1 


0.  r> 


Fig.  6:  Comparison  of  RMS  derived  from  sample  size  10. 

(Step  function  estimators  vs.  Bayesian  estimators) 


56 


. fa  1 H ..  l.M.  -• 


*>,'M.i  .V‘.*  PMH  .vs.  .■‘><v».  ,vs,  t*M\) 


I ..»  .Hh  - Jl 


>- » i2i  - ): 


».  \ - ifc 


*.  1 is*  - V 


/.>Sit - }Z 


»•  ^ .'»L  - n 


j. ) 


J.O 


•’•  ’ C.  f'} 


-M 


■: . ( ( ‘r -0-' 


f . ; rcf  _,- j 


f . 1 1 '(--O’ 


?.  ’<il'  -J2 


( i ft  --;> 


I .( i; 


Fig.  7|  Comparison  of  rms  derived  from  Scunple  size  25. 
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Fig.  8;  Comparison  of  RMS  derived  from  sample  size  25, 

(Step  function  estimators  vs.  jackknife  estimators  of 
logistic  transformation) 
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APPENDIX  A 


ESTIMATORS  FOR  GROUPED  DATA 


With  grouped  censored  data  the  definition  of  p(t^/t^  given 
by  equation  (5)  does  not  hold  unless  the  assumption  is  made  that  all 
truncations  occur  at  the  end  of  the  time  interval.  If,  on  the  other 
hand,  it  is  assumed  that  all  truncations  occur  at  the  beginning  of 
At^  the  equivalent  form  of  equation  (5)  is 

N.  - a.  - r. 

P(t./ti_i)  = (G-i) 

1 1 

Witli  elements  were  presents  at  beginning  of  interval,  i.e.,  at 
time  r^  elments  failed  during  the  interval,  and  a^  elements 

truncated  from  the  sample  during  the  interval  but  prior  to  failing. 

As  a hypothesis,  assume  that  all  aborts  occur  simultaneously  somewhere 
within  the  time  interval,  so  tliat  r'  failures  occur  prior  to  the 
truncations  and  time  remaining  r^  - r ' after  the  truncations.  Then 


P(t./ti_i)  = 


N.  -r'  N.  -a.  -r. 
1 ^ 111 

N.  ’ N.  - r.  - a. 

1 111 


Thus,  the  value  of  P(t^/t^_^)  depends  on  when  the  truncations  occur. 

It  is  assumed  that  this  is  not  known  for  the  grouped  data  case.  Never- 
theless, it  is  possible  to  place  limits  on  the  value  of  p(t^/t. 
since  equation  (G-2)  always  gives  values  between  those  of  equation  (5) 
and  (G-1) . Thus 


N,  - a.  - r.  N.  - r. 

— < p(t  /t  ) < — i- 

N.  - a.  _ V i-l'  - N. 
11  1 


Tor  average  sample  size  approximation,  a simpler  expression  from  the 
point  of  view  of  computational  ease  may  be  derived  by  substituting 
a/2  for  a in  equation  (G-l)  giving 
tl  - I - r 

p(t.  ) = (G-4) 

X - a 

j - j 

The  equation  (G-4)  may  be  thought  of  as  the  result  of  assuming  that 
the  average  number  of  elements  in  the  time  interval  is  the  number  at 
the  beginning  decreased  by  half  the  number  of  truncations. 

Records  are  usually  available  to  provide  a fairly  precise  time  the 
deatii  events.  In  the  medical  example,  the  exact  time  of  death  is 
usually  recorded  in  medical  records  required  by  law.  In  the  equipment 
lifetesting  example,  the  time  of  malfunction  or  failure  is  usually 
known  very  precisely  if  the  results  are  catastrophic;  and  maintenance 
records  give  a reasonably  precise  time  even  if  the  failure  is  not 
critical  to  a larger  system.  In  the  military  example,  the  event  of 
interest  is  usually  a sensor  detection  or  some  other  action  that  is 
routinely  recorded  in  a log  book. 

Equaiton  (G-4)  is  a modification  to  the  product-lirait  estimator, 

P2»  when  the  times  of  truncation  are  known  only  in  grouped  form.  Herd, 
reference  (2),  suggests  a similar  modification  to  estimators  using  the 
second  approach  (p^  or  p^)  with  aggregated  truncation  data.  Illustrate 
results  for  this  method  based  on  the  sample  data  base  of  the  main  test 
are  given  below.  Here,  of  course,  we  do  not  know  that  individual  B 
dropped  out  of  observation  at  time  2 and  that  individual  D dropped  out 
at  time  6.  We  know  only  that  the  two  truncations  occurred  in  the 
interval  (1,3)  and  (3,7),  respectively. 
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Product  limit's  modification  is  denoted  by  p^'Ct)  and  Herd's  modi- 
fication is  denoted  by  p^. ' (t)  • 

Their  results  on  the  sample  data  base  are  as  follows. 


t 

0-1 

1-3 

3-7 

(7) 

t 

0 

1 

3 

7 


P2*(t) 


5/5  = 1.0 
4/5  X 1.0  = 0.8 

2. 5/3. 5 X 0.8  = 0.571 

0.5/1. 5 X 0.571  = 0.190 


P3' (t) 


1.0 

5/6  X 1.0  = 0.833 
3. 5/4. 5 X 0.833  = 0.640 
1.5/2. 5 X 0.648  = 0.389 


APPENDIX  B 


LISTING  OF  COMPUTER  PROGRAM 

C TRUNCATED  C/STA  PRCGRAM 

C 

01  MENS  IGN  NN(9)  , D(  14)  ,T(  100  j),  IT  ( KOO)  ,Pl(lC:C),Pi(lCC 
=*0  ) ,P3(1  000)  ,P4(10C0)  , P5(100O),P6(x  JOO),  PJA4aC0j)  ,PJA5 
( 1000)  f PJA6  11000  ) ,TJ  ( 1000)  ,1TJ(1000),P2(10JC),S(14,9), 
^0  (14,9)  ,W(  14,9),  PHI  1000)  fPli(10CC),P13(  lOOC  ) , SL 1 < 50 , 5 
*C) , SLZ( 50 ,50 ,SL5  ( 50,50) ,S3AR1 (30  ) ,S3AR2 (50  ),SEAS3 (50  ) 
>i-P 04(50, 50)  ,PJ51  50,50)  ,P06(  50,50)  ,PJ2( 5 0,50)  ,F22(  1 COO  , 
) ^SL4(50,50)  ,SBAR4(  50)  , PJA2(1  jjO  ) , C ( I't , 9 ) , RM  S 1 ( 50 ) , RMS2  ( 

=!'(  50)  ,RMS5(  50)  ,RMS4(50)  ,UINT1  (50)  ,LINT2  (50)  ,LINT^(50),U 
*INT3(50),RINTl(5J),RINT2(50),RlNT5(50),RINT4(50),CF(e) 
*VARJ2 (50) , VARJ4( 50) , VARJ5(50 ) ,VARJ6(50),UPi(5C),UP2(5C 
^)  ,UP3(50) ,UP4( 50) ,ROl{50) ,Ra2(50)  ,R03(50)  ,RC4(50) 

CALL  CVFLOW 
C 

1 FORMAT  (I5,I10,F10.4,5I5,F3.3) 

2 FCRMAT( IX, • NOTE  ERROR') 

3 FCRMATdX,' NTRUNC  ERROR') 

4 FCRMAT( IX, 'NREPL  ERROR') 

5 FORMAT ( IX, • NEVENT  ERROR') 
t FCRMAT( 1X,'NWRITE  ERROR') 

7 FCRMAT  (3X,9F6.3) 

H FCRMAT  ( 1X,2I5 ,l4Fe.5  ) 

9 FCRMAT  (IX, 'P'  ,I 2,9Fe.3,3X, 'MEAN'  ) 

1C  FCRMAT ( ' • ) 

11  FCRMAT( •!'  ) 

12  FCPMAT( IX, I 5,10F1C.5) 

13  FCRMAT  (IX,  2 1 5 , F 10 . 4, 4 1 5 ) 

14  FORMAT  ( 1 X, 'P'  , I 2 ,9F8.3 ,3X,  ' ABS' ) 

15  FORMAT  (6X,9(14,4X)  ) 

16  FCRMAT  (IX, 'P'  ,I2,9F8.3,3X,  'RMS' ) 
la  FCRMAT  (IX, 215) 

C 

C READ  INPUTS  AND  SET  INITIAL  VALUES 

C 

A=C.01 

2 5 READ(5, 1 ) ND IE , NT RUNC , X TkUNC , NR E F L , NEV ENT , NW F I T E , ISEEC, 
=*'NENO,TN 
JEVENT=NEVENT-1 

V^RITE(o,13)NDIF  ,NTRUNC,XTRUNC,NREFl,NEVENT  ,NWR  ITE,  ISEE 
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IF  (NDIE.LT.l)  GOTO  105 
IF  INDIE. GT. 2)  GOTO  105 
IF  ( NTRUNC. LT.l ) GC  TO  105 
IF  (NTRUNC. GT. 3)  GOTO  103 
IF  (NREPL. LT.l ) GO  TO  104 
IF  (NREPL.6T.1000)GOTO  104 
1F(NEVENT.LT.2 ) GOTO  102 
IF  (NEVSNT.LE.1000)  GCTQ  200 
C 

C EFRCR  MESSAGES 

C 

102  WRITE  (6,5) 


C 

C 


WRITE  (6,10) 

NP  = 14 

SRH  = SURT( .5  ) 

DC  26  1=1,1000 
T(I )=0. C 
Pl(  I)=0.0 
P2(  I)=0.0 
P3(  I ) = 0.0 
P4(  I ) = 0 .0 
P5 ( I )=0.0 
Pfcd )=0.0 
TJ(  I)=0.0 
PxK  I )=0.0 
P12(I )=0.0 
P13(  I)  = 0.0 
26  P2( I )=0.0 

TEST  INPUTS 
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103 

STCP 

WRITE 

( t,3) 

iC4 

STCP 

WRITE 

( 6,4) 

105 

STCP 

WRITE 

(6,2) 

ICE 

STCP 

START 

MAIN 

2CC 

CC  250 

J = I, 

NN( J) = 

0 

CC  250  1=1, NP 
S(I,J)=0 
U(  I , JJ  =0 
C(  I,J)  = 0,0 
25C  I,J)=3 

DC  4999  IREFL=1 ,N9EPL 
NCI=0 

DC  999  ieVEM=l  ,NEVCNT 

CREATE  TTRUNC( ) 

CALL  RANDOM! ISEED,TTR,1) 

GCTQ  (300,350,400 JtNTRUNC 
GCTG  103 

5CC  TTR=TTR*XTRLNC 
GCTG  500 

3fC  TTR=-XTRUNC*ALOG(  TTR ) 

GCTG  500 

4C0  TTR  =XTPLNC 

CREATE  TDIE(J 

5CC  CALL  RANDOM!  ISEED,T0I,1) 

GCTG  (800,700) ,NDIE 
GCTG  102 

700  TCI=-AL0G! TCI ) 

DETERMINE  SMALLER  CF  TDIEO  ANC  TTRUNC!) 

8CC  IF  (TDI.LE.TTR)  GGTO  310 
TT=TTR 
ITT=0 
GCTO  850 
81C  TT=TOI 

NCI=NDI +1 
ITT=1 

aSC  T!  IEVENT)=TT 
IT!  lEVHNT)  = ITT 

CRCER  DATA 

67C  I F (lEVENT.EQ.l ) 

II=1EVENT-1 
CC  890  1 = 1 , II 
IF  (TT.GT.T!!)} 

I 11=11-1  + 1 
CC  630  J = 1 , III 
JJ=IEV£NT-J 
I T( JJ+1 J=IT!JJ) 

3eC  T(JJ+1 )=T( JJ) 

IT! I )=I TT 
T ( I ) = TT 
GCTG  999 
aSC  CCNTINUE 


GCTG  999 
GOTO  890 


999  CCNTINUE 

TMsT!  N’EVENT  ) 

T7  = 0 

DII=NDI+2. 

IF  (NOI  .GT.O)  GOTO  8122 
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P^(NEVENTJ=SRH 
PrCNEVENT )=SPH 
TT=0 

DC  8121  J=1,N£VENT 
8121  TT=TT+T(J) 

DN=TT/r (^EVENT ) 

P£(NEVcNT)=SQRT(DN/(CN+l.J ) 

GCTC  1111 

CALCULATE  Pl()  AND  P4()  AND  PIK)  VECTORS  AND  P7  DATA 


I 11=0 

DC  2199  1=1,NEVEMT 
T7  = T7+T  ( I ) 

IF( IT( I ).Eg.l)  GCTC  2150 
J = 1 

GCTC  2199 

21 5C  Pl( I)=FL0AT(N-1 ) /float (NDI) 

PA( I )=PL0AT(N)/FL0AT(NDI+1) 

PIK  I )=FLCAT(N)/DI  I 
TTT=T( I ) 

111=11 
I I = I 
J = C 
N=N-1 

2199  CONTINUE 

T 7=-NDI /T7 

IF  IJ.EC.O)  GOTO  2221 
TI I=T( I I) 

CTI=T(I  in-TII 
IF  (NDI.GT.l)  GOTC2200 
TAL=-TI I/ALCG(  PA( I I)  ) 

TTAU=-TII/AL0GIP11  (ID) 

GCTC  2210 

22CC  TAL=  DTI/ALCG( PA( I I ) /PA( 1 1 I ) ) 

TTAU=(  T(  I I I J-TI  I )/ALaG(PlK  1 1 )/Pll  ( T 1 1)  J 
221C  DT=TN-TII 

IF  (OT.GT.ISO’I'TAU)  TAU  = 0T/i50 
I F(DT  .GT.150^TTAU)TTAU  = iJT/150. 

Pli(NEVENT) =Pii  ( I I )^EXP(-OT/TTAU) 
PA(NEVENT)=P4(II )* EXP ( -DT/T AU) 

2221  IF  (NCI.NE.NEVENT  ) C-0  TO  2225 
DC  2224  I=1,NEVENT 
F2( I )=P1( I ) 

P3(  n=PK  1 ) 

P12(  I)=P11(  1) 

Fl3( I)  = P11(  I) 

F£(I)=P4( I ) 

2224  PE( I ) = P4( I ) 

GCTC  nil 

CALCULATE  P2()  AND  P5()  AND  P12()  VECTORS 


2225  PF=1. 

N=NEVENT 

PFF=FLaAT(N+l)/FLCAT(N+2) 

P = l. 

J = C 

DC  2399  I=1,NEVENT 

IF  ( IT( I ) .EQ.l ) GCTC  2350 

J = 1 

GCTC  2399 

235C  PP=PP*FL0AT(N-1)/FL0AT(n) 
P=P*FLOAT(N )/FLOAT(N+l J 
PFF  = PPP«PLCAT(N)/FLGAT (N-n  ) 
P2( I )=PP 
F5(  I )=P 
P12(D=PPP 
J = C 
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23=9  N=N-1 

IF  (J.cfi.O)  G0TC2425 
IF  (NOI.GT.l)  GOTO  2400 
T3L=-TI I/ALCG( P5 ( II  ) ) 

TTAU=-T I I/ALCG( P12  (I  I ) ) 

Gi;TC2410 

24 CC  TAL  = 0TI /AL0G(P5 ( I I )/P5  ( 1 1 I ) ) 

TTAU=DTI/ALCG(P12{ 1 1 ) /P12( I 1 1 ) ) 

241C  IF  (DT.GT.150*TAU)  TAU=DT/150 
P5(NeVtNT)  =F=*EXP(-0T/TAU) 

I F(0T.GT.15C*TTAU  >TTAU=OT/150. 

P12(NEVENT) =FPP*EXF(-CT/TTAU  ) 

CALCULATE  P5()  ANC  P6()  AND  P13()  VECTORS 

2425  PP=1. 

N = NEVEiMT 
P = l. 

PFP=FLOAT{iN  + l)  / float  (N+2J 

J = C 

TT=0 

TTT=0 

DC  2599  I=1,NEVENT 

IF  ( IT(  I ) GOTO  2550 

J = J + 1 

TT=TT+T ( I J-TTT 
GOTO  2599 
255C  0N=N 

IF  (TT.NE.O)  0N=0N  + TT/(T( I 1-TTT  ) 

PF=ppc( oN-1 )/ON 
P=F*DN/ (DN+1) 

FFP  = PPP*CN/ (ON-H  . ) 

P2(I )=PP 
P6 { I )=P 
P15( I)=PPP 
J = C 
TT  = 0 

TTT=T( I ) 

2599  N=N-1 

IF(J.E0.0)GG  TO  llil 
CTT=TN-TTT 

IF(CTT.GE.1E-70)GCTC  3005 
CA=.5=!‘(  J + 1 ) 

GOTO  3010 
50C5  CN=TT/OTT 

30  1C  P£(NEVENT)=P*SGRT(CN/ (DN+1. » ) 

IF(\0I.GT.l  )GO  TO  1997 
TTAU=-TII/AL0G{P13(II  )) 

GOTO  1998 

199?  TTAL=DTI/ALCG(P13(II)/P13(I  ID  J 
1998  IF ( OT. G T. 1 5C*TTAU) TTAU=DT/1 50. 

P13(NEVENT )=PPP*EXP(-CT/TTAU) 

SET  LP  A LOOP  FOR  ALL  JACKKNIFE  CALLS  (P J4, P J5 , FJ6 ) 

1111  CC  1000  1=1 tNEVENT 
DC  1011  J=1,NEVENT 
PJ2(I, J»=O.C 
PJ4( I , J)=O.C 
Pj5(I, J)=0.C 
PJ6(I,J»=0.C 
SL1(I,J)=0.C 
SL2(I»JJ=0.C 
SL3( I , J)=0.C 
SL4( I, J)=O.C 
1011  CCNTINOE 

PJA2( I ) =0.0 
PJA4(I )=0.0 
P-A5(I)=0.0 
PJA6( I ) =0.0 
SEARl ( I i=O.C 
SEARZl I J=O.C 
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S£^R3( I )=0.C 

SeAR4( I ) =0  .c 

P Z 2 ( I)  = 0 . J 
VARJ2( I )=0.C 
VARJ4( I )=0.C 
VARJ5( I }=0. C 
VARJ6( 1 )=0.C 
lOCC  CONTINUE 

DC  1001  I=l,NeVENT 

MOVES  DATA  INTC  TJ()  AND  ITJd 

K = 1 

JNEXT=0 

jeEFCR=0 

JAFTER'O 

10C2  IFIK.EU.I )GC  TO  7003 
TJ(K)=T (K) 

I IJ(K)  = IT(K) 

IF(  IT(K  I .EQ.OJGO  TC  7001 
JNEXT=JBEFCP 
jeEFCR=K 
70C1  K=K+1 

GC  TO  1002 
70C2  JAFTER=J8EFCR 
jeEFCR= JNEXT 
GC  TO  1010 

70C3  IF(I.GT.JEVENT)  GO  TC  7002 
1003  IFIK.GT.JEVENT JGO  TO  1004 
TJ(K)=T (K+1 J 
nj(K)  = lT(  K+1) 
IF(ITJ(K).EC.O)GO  TO  4002 
IF( JAFTER.EC.O ) JAFTER=K 
40C2  K=K+1 

GC  TO  1003 

ICC  A IFl JAFTER  .EC.O) J AFTER  = J EVENT 
lOlC  NOIJ=NDI-IT(I) 


CHECK  IF  ZERO  DEATHS 

I F (NCIJ  .EG.C  )GOTC  1001 

N=NDI J 

P=l. 

J = C 


I I = C 
I 11  = 0 

GC  TC  1014 

CALCLLATE  PJ4(  )VECTORS 

1014  DC  lOlo  IJ=1,JEVENT 

I F ( ITJ(  IJ) .EQ. 1 )GG  TO  1015 
J = 1 

GC  TO  1016 

1015  P=FLOAT(N )/FL0AT (NCIJ+1  ) 

PZ( I J) =P 

111=11 
I I = IJ 
J = 0 
N=N-1 

1016  CCNTINUE 
lF(J.Eg.0)GC  TO  1019 
T 1 IsTJl  1 1 J 

CTI  = TJ(  Iin-TI  I 
IFINDIJ.GT.DGO  TC  1017 
TAU=-TI  I/ALCG(  PZ(  ID  ) 

GC  TO  1018 

1017  TAL  = DTI /AL0G(P  Z( I I )/PZ(  II  IJ  ) 
IJie  DT=TJ( JEVENT)-TI I 

IF(DT.GT.150.#TAU )TAU=DT/150 
PZ(  JE VENT)  =F’!'EXP  (-DT/ TAJ) 
1019  K=1 


VECTORS 
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iOiC  IFCK.EU.IJGC  to  2021 

PJ^(K,  I ) = ALCG{(  PZ(K)+A)/(1-PZIK  )+/S)  ) 

K = K + 1 

GC  TO  1C20 

2021  IF(K.EQ.NEVENT)GO  TO  5021 
IF(IT(I ).EQ.0)G0  TC  1025 
5021  TJA=TJ( JAFTER) 

IFUBEFCR.NE.O)GO  TO  1022 
PX  = ALOG( FZ( JAFTER  ) )*T ( IJ/T JA 
PZZ=0.0 

IF  (PX.LT.-150)  GOTO  7025 
PZZ=EXP(PX) 

7025  PJA(I,I)=ALCG((PZZ+A)/(1-PZZ+A) ) 

GC  TO  1325 
1022  7«B=TJ( JBEFCR) 

DT J=TJA-TJ8 
CTE=T(I)-TJE 

PX  = DTB#ALOG(PZ( JAFTER )/PZ( J6EFCR  ) )/DTj 
PZZ=0.0 

IF  (PX.LT.-150)  GOTO  7026 
PZZ  = PZ( JBEFCR)*EXP(PX  ) 

70 2 1 PJA(I,I ) = ALCG( (PZZ  + A)/(l-PZZ  + Ai  ) 

1025  IFIK.GT .JEVENTIGCTC  2C27 

PJA(K  + 1 ,I }=ALOG(  ( FZ(K)+A)/(1-PZ(K)+AJ ) 
K = K + 1 

GC  TO  1025 

3027  IF(NDI.NE.NEVENT)GGTG  1026 
DC  3028  K=1,NEVENT 
CC  3028  L=i,NEVENT 
FJ5(K,L)=PJA(K,L) 

5028  PJc  (K.,L  ) = PJA(K,  L ) 

GOTO  iOOl 

CALCLLATE  PJ5 ( ) VECTORS 


1026  N=JEVENT 
P=l. 

PP  = 1. 

J = C 

CC  1028  IJ=1,JEVENT 

IF{  ITJ<  IJI.EQ.DGC  TC  1027 

J = 1 

GC  TO  iC28 

1027  F = P=»FLOAT(N)/FLOA-^(N  + 1) 

FF  = PP=«=FLCAT  (N-1  )/FLO/!T(N) 

PZ( I J)=P 

PZ2( I J )=PP 
J=0 

1026  N=N-1 

IF(J.EO.O)GC  TO  1C31 
1 F(NDIJ  .GT.l  )G0  TC  102v 
TAL=-TI I/ALCG( PZ ( 1 1 ) ) 

GC  TO  1330 

102  9 TAL=DTI/ALCG(PZ( I I )/FZ(  II  1)  i 
10 'C  IF(DT.GT.150*TAU) TAU=DT/150. 

PZ( JEVENT )=F*EXP{-CT/TAU) 

1031  K=1 

1032  IF(K.EU.I )GC  TO  2033 

PJ5(K,,I  ) = ALCG((PZ(K>+A)/{  1-PZ(K  » + A)) 
PJ2(  Kf  I )=ALCG( (PZ2(K)+A)/(1-PZ2(K)+A)) 
K = K + 1 

GC  TO  1032 

2033  IF(K.Eg.N£VENT)GG  TO  5053 
IF( IT( I ).EQ.O)GO  TC  1136 
5033  IF( JBEFCR.NE.O JGO  TO  1135 

PX  = T(  DIALOG  (P2(  JAFTER)  J/TJ  A 
PZZ=0.0 

IF  (PX.LT. -150)  GCTO  7027 
PZZ^EXP (PX) 

7027  PJ5(I,I )=ALCG( ( PZ Z+ A ) / ( 1-PZ Z ♦A ) ) 

GC  TO  1136 


1135  PX=DTB*ALOG(PZ( JAFTER )/PZ(J6EFCR) )/DTJ 


PZZ=0.0 

IF  (PX.LT.-150)  C-CTO  7028 
FZZ=PZ( jeEFCR)*EXP(PX  ) 

70 2 £ PJ5( I ,I  )=ALCG( ( PZZ  + A) / ( 1-PZZ  + AI  ) 

1136  IFU.GT.JEVENT  )G0  TO  1036 

PJ5(K+1,I)=ALCG((PZ(K)+A)/(1-PZIK)4A)J 
P J2(K.  + 1 , 1 >=ALOG(  (PZ2(Ki+A)/(l-P22{K)4-A)) 
K = K + 1 

GC  TO  1136 
C 

C CALCLLATE  PJ6()  VECTORS 
C 

103c  N=J£VENT 


T7  = C 
TTT  = 0 

DC  1033  IJ=1,JEVENT 

IF(  ITJ(  IJ) .EQ.l )GC  TO  1037 

J = 1 

71=TT  + T J( IJ  )-TTT 
GC  TO  1058 
1037  CN=N 

IF(TT.NE.O) CN=DN  + TT/ (TJ ( IJ ) -TTT  ) 

P = P=t‘DN/  (0N>1.) 
pz  n J J = P 
J=C 
T7  = C 

77T=TJ(  IJ J 
1036  N=N-1 

IF(J.EQ.O)GC  TO  ICAl 
CT=TJ( JEVENT)-TTT 
IF(DT.GT.1E-70)GC  TO  1039 
0N=.5*( J+1) 

GC  TO  1040 
1C3C  DN=TT/07 

IC^C  PZ(JEVENT)=F=!<SQRT(CN/(  DfNJ+l.  ) ) 
lOAl  K=1 

10^2  IF(K.EQ.I)GC  TO  2043 

PJ6(K,I)=ALCG((PZ(K)+A)/(1-PZ(K)+A)) 

K=K  + 1 

GC  TO  1042 

2043  IF(K.EO.NEVENT)GO  TO  5043 
IF(  IT(1 ).Ee.O)GO  TC  1146 
5042  IF( J6EFCR.NE.0)GC  TO  I0t5 

PX  = T( I ) 4AL0G(PZ( J AFTER)  )/TJA 
PZZ=0.0 

IF  (PX.LT.-150)  GOTO  7029 
FZZ=EXP(PX) 

7C2S  Pj6(I  ,I  »=ALCGH  PZ  Z+A  ) / ( i-PZ  Z +A ) ) 

GC  TO  1146 

1045  PX=DTB*AL0G(PZ(  J AFTER  )/PZ(  J8EFCR  ) )/'DTJ 
PZZ=0.0 

IF  (PX.LT.-150)  GOTO  7030 
F2Z  = PZ( JBEFCR)*EXP  (PX  ) 

70 2C  PJ6( I ,I )=ALCG{ ( P Z Z+ A ) / ( 1-PZ Z+A) ) 

1146  I F(K.GT.jEVENT)GOTO  1001 

P J6(K+l ,I )=AL0G( ( FZ(K )+A)/ ll-PZ IK )*A) ) 

K = K + 1 

GC  TO  1146 
lOCl  CONTINUE 

DC  8888  I=1,NEVENT 
DC  8811  J=1,NEVENT 
I F( I .EQ.J)GCTO  8811 

SLl(IfJ)=FLCAT(NEVENT)*AL0G((P4(I)+A)/(l-F4(I)+A))-FLC 
*AT(NEVENT-1  )*PJ4l  I,J» 

SL2(I»J)=FLCAT(NEVENT)*AL0G((P5(I)+A}/(i-P5(I)+A))-FLC 
*AT (NEVENT-l  )*PJ5(  I , J ) 

SL3( I , J )=FLCAT (NE VENT )#ALQG( (P6(I)+A)/(l-Pe(!)*A})-FLC 
*AT(  NEVENT-1  )>^PJ6  ( I , J) 

SL4( I,  J)  = FLCAT(  MEVENT )«AL0G( <P2( I ) +4 ) / ( 1- P 2 ( I ) + A » ) - FLC 
■*AT(NEVENT-i  )’f‘PJ2  ( I , J) 
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d8il  CCMINJE 
SacS  CCNTINLE 

DC  9999  I=1,NEVEM 
DC  9911  J=1,NEVENT 

Sa^Rl ( I )=SL1( I,J)/NEVEMT+SdAR1 (I  ) 

SeAR2(  I i = SL2( I , J) /NEVeNT+S6AR2( I) 

SeAR3 ( I ) = SL*( I, J )/NEVENT+SBAR3( I ) 

Se  AR4(  I ) = SL4(  I , J ) /NE  VENT<-SBAR4  ( I ) 

9911  CCNTINUE 

IF(SBAR1(  ! ) .LT.-18D.  ) SBARK  I)=-1£C. 

IF(SBAR1(  I)  .GT.17A.)  SBARKI  )=17A. 
iF(S8AR2(  n .LT.-iec.  )SBAR2(  n=-iec. 

I F (S8AR2(  1)  .GT .1 7A . ) SBARE ( I ) = 174 . 

I F(SBAR  2(  1 ) .LT,-1 tC. ) SBAR5( 11=- ISC. 

IF(SBAR3( I) .GT .174. )SeAR3( I }=174. 

IF<iBAR4( 1)  .LT.-l 80. )SBAR4( Ii=-160. 

I F( S8AR4( I ) .GT .1 74. ) SBAR4( I 1=174. 

EL1=EXP(SBAR1 (111 
EL2  = EXP( S3AR2( I 1 1 
EL3  = EXP (SBAR3( I 1 ) 

EL4=EXP( SBAP4( I ) 1 
PJA4( 1 1 =EL1/{1 . + EL11 
PjA5 (I  1 =EL2/(1.+EL2 ) 

PJA6  (I  1 =EL3/(1 . + EL3) 

PJA2  (I  1 =EL4/(1.+EL41 
999?  CCNTINUE 

DC  8999  I=1,NEVFNT 
DC  8911  J=1,NEVENT 
I F { SL4(  I , J1  .EQ.C  .0 IGOTO  3111 

VARJ2(I 1 = (SL4(  I , J 1-SBAR4( I 1 1**2/FLCAT(NEVENT-11+VARJ2( 
1 

3111  IF(SL1 ( I ,J1  .EQ.O.OIGQTO  3112 

VARJ4(Il  = (SLl(I,J)-S3ARi(I}  1 ##2 / F L C AT ( NE V EN T- 1 ) + V A R J4 ( 
*11 

3112  IF(SL2( I, Jl  .EO.O.OIGOIO  3113 

VARJ5(I1  = (SL2(I,J)-SBARZ(I11=)‘*2/FLCAT(N£VENT-1)+VARJ5( 
^ I J 

Sll'^IFCSLSt  I , J)  .EQ.O.OGOTO  8911 

VARJ6  II 1 = ( SL3(  I , J 1-S3AR3( 1 1 1**2/FLCAT(NEVENT-1 1+VAR je ( 
*11 

89*1  CCNTINUE 

RNS1(I)=SQRT(VARJ2(I ) /F lOAT ( NE V EN  1 ) 1 
RRS2 ( I 1 =SURT(VARJ4( I ) /FLOAT  I NEVENT ) 1 
RRS3(I)=S.jRT(VAPJ5(Il  /FLOAT  ( NE  VENT  1 ) 

RyS4(I  1=SQRT(VARJ6(I  l/FLQAT(NEVENTn 
UPl(Il=SBAR4(Il+PNSim*TN 
R01(I)=SbAR4(I)-RRSl(Il*TN 
UF2(I1=SBAR1(I)+RMS2(I)*TN 
RC2( I 1 =SB ARi { 1 1-RRS2 ( I )*TN 
UF3(I)=SbAR2(Il+RNS3(I 1 *TN 
RC5(I)=SBAR2(I)-RyS3(I)*TN 
UF4( 1 1 =£BAR3 ( I ) +RyS4( I ) *TN 
RC4(I 1=S8AR;( I)-RyS4{ I 1*TN 
IFIUPK  I ) .GT.174.  lUPK  I 1 = 174. 

IF(UPi( I1.LT.-100.1UP1(I1=-100. 

IFIROK  I1.GT.174.)RG1(I)=174. 

IFIROK  II  .LT.-IOO.  IROK  II  =-100. 

IF (UP2( I 1.GT.174. 1UP2( 1 1 = 174. 

IF(UP2( I). LT.-IOO. )UP2( 11=-100. 

IF(R02(  I ) .GT.174. )R02 ( I 1 = 174. 

IF(R0  2(  1 1 .L  T.-IOO.  1R02( I 1=-100. 

IF (UP5(  I ) .GT.174. 1UP3 ( I 1 = 174. 

IF(UP3(  I1.LT.-10C.1UP3(I1=-100. 

IF  (R03(  1) ,GT.l 74. 1R03 ( 1 1=174. 

IF (R03(  1 1 .lT.-IOO. 1R03 ( 1 1=-100. 

IF(UP4{  n.GT.174.1LP4(Il  = 17  4. 

IF(UP4< I).LT.-100.)UP4(I)=-100. 

IF(RO^(  I I .G 7.1 74. 1R04( 11=17  4. 

IF(R04(  I 1 ,LT.-100.  )R04( Il=-100. 
LINT1(I)=1./(1.+1./EXP(UP1(I11} 
RINTl(Il=l./(l.+l./FXP(Rai(Il)l 
UI NT 2(11  = 1. /(I. +1./EXP( UP  2(111) 
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RINTid  >=1./(1.  + 1./EXP(R02(I))) 

UIM3(  I)  = l./(l.+l/EXP(UP3(in) 

RIM5(n  = l./(i.+l./EXP(R33(I))) 

UINT4(I J=1./(1.+1./EXP(UP4(I))) 

RINT4( I )=1./(1 .+1  ./EXP(R04(  1) ) ) 

89EE  CCNTIN'JE 

PRINT  OUTPUT 

3050  IF  (NKRITE.EQ.OJ  GCTC  3500 

WFITE(6,S)(IREPL,IT(n,T(n,Pim,F2(n,P5(I)»P4tn.P3 
♦(  1)  tP6(  I)  , PJA2(  I ) ,PJ/S4(  n , PJA5  ( I)  »FJA6(  I ) ,P11  ( H , Flh  I 
*)  ,P13( I ),I  = 1,NEVENT) 

WRITE  (EtlOJ 

RECLCE  VECTORS 

35JC  K=1 

CC  40&0  I=\,NEVENT 
IF  ( IT{ I I.EC.l)  GCTC  5900 
IF  ( I .NE.NEVENT ) GCTO  AOJO 
B9CC  T (KJ=T(  I ) 

Pj.(K)=Pll  I ) 

P2(K)  = P2(  I ) 

F3(K)=P5(  I ) 

P4(K) =P4( I ) 

P5(K)=P5(I) 

P6 (K)=P6(  I ) 

FJA2<  K) =PJA2(I ) 

FuA4( K) =PJ A4( I ) 

PwA5 (K) =PJA5( I ) 

Fj46(  K) =PJA6  (I  ) 

P11(K)  = P11(  I) 

P12(K) =P12(  I) 

P13(K)=P13(  I ) 

UINTKK  )=UINT1  ( I ) 

RINT1(K,)=RINT1  I I ) 

LINT2(K )=UINT2( I ) 

RINT2tR>  = RlNT2d  ) 

UINT3(k) =UI NT3 ( ! ) 

RINT3(K )=RINT3(  I ) 

UINT4(K )=UINT4(  I ) 
kINT4(R)=RINT4( I ) 

K = K + 1 

4CCC  CONTINUE 


C/!LCULATE  DIFFERENCES,  i4EAN,  RMS,  AND  MEAN  AbS 


K = 1 
TTT=0 
CT=-T(l ) 

PP1  = 1 
PF  2 = 1 
FF3=1 
PF4=1 
FF5  = 1 
PF6=L 
PF6  = 1 
PF9=1 
FF1C=1 

PFll=(NCUl.)/niI 

PF12= (NEVENT+1.  )/ (NEVENT  + 2.  ) 

PFi3=PP12 

FFI4=1. 

F;4=P4( 1) 

P(,5=P5(  1) 

FCe=P6( I ) 

FC0=PJA4(1 ) 

PC9=PJA5( 1) 

PC10=PJA6 (1  ) 

CCNF1=1  . 

CCNF2=1  . 
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CCNF3=1  . 

CCNF4=1  . 

CCNF5  = 1 . 

CCNF6  = 1 . 

CCNF7=1  . 

CCNF8-1  . 

T4=DT/4L0G(Ft4) 

t;=dt/alog( PCS ) 

Tfc=0T/ALCG(FC6) 
ie=DT/ALOG{ FC81 
Tc=DT/ALOG( PQ9) 

T1C=DT/ ALCG (PQIC ) 

DC  4998  J=lf9 

TT  = .1>*'J 

P=1-TT 

IF  (NDIE.EQ.2)  TT=-ALaG(P) 

IF  (TT.GT.TN)  GOTC  4999 
NM  J)=NN(  J)  +1 

4iCC  IF  ITT.LE.KK})  GCTC  4200 
PF1=P1 ( K) 

FF2=P2( K) 

FF2=P2( K) 

FF4=PQ4 

PP5=Pw5 

FF6=PQ6 

FF8=PQ3 

PF9=PC9 

PF1C=PQ10 

PFll=Pil(K) 

PF12=P12( K) 

PFIO^PI 3(K ) 

PF14=PJA2 (K ) 

CCNF1=UINT1 (K) 

CCNF2=RINT1 (K) 

CCNF3=UINT2 (K) 

CCNF4=RINT2(K) 

CCNF5=JINT3 (K) 

CCNF6=R INT3 (K) 

CCNF7=U INT4(KJ 
CC\F8=R  INT4 (K) 
nT  = T(K) 

K = K + 1 
PC4=P4( K) 

PC5=P5( K) 

PC6=P6( K) 

PC8  = PJA4<K  ) 

PC9=P JA5(K) 

PtlO=PJA6(K) 

DT=T(K) -TTT 
TT1=ALCG( PPE/PC8) 

TT2=AL0G( PF9/PQ9 ) 

TT2=ALOG( PFIO/PQI C) 

IF(  P08.EQ.0.C.0R.  TTl  .EQ.0.0  J TT1  = ALCG(  I PP8+A)  / ( 1-P(.8  + Ai) 

IF  (PC9.  EU.0.0.0R.TT2.EQ,0.0  )TT2=ALCG(  ( PP9+A  ) / ( 1-PU9  + A )) 

IF  (PgiO.EQ.C.O.OR.TT3.EQ.O.O)TT3  = ALOG(  ( PPIO  + A)  / (1-F(.1C  +A)) 
T4=0T/ALCG( PP4/PQ4) 

T5=DT/ALCG( FF5/PQ5 ) 

T6=DT/ALCG( FP6/PQ6) 

Te=DT/TTl 
T9=0T/TT2 
T1C=DT/TT3 
GCTC  4100 
42CC  0(1)=FP1-P 
C(2)=PP2-P 
C (3  ) = FP3-P 
C17=TTT-TT 

C(4)=PP4*EXF(DTT/T4)-P 
C (5)=PP5#EXF(0TT/T5)-P 
D( 6)=PP6*EXF(DTT/T6)-P 
D(7)=EXP(TT#T7)-P 
C(6}=PP8*EXF(OTT/T8)-P 
D( 9)=PP9#EXF(DTT/T9)-D 
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4911 


4997 

499£ 


C (10)=PP10=«‘EXP(  DTT/T10)-P 

D(  11  )=PP11-F 

Ca2)=PPl2-F 

C(13)=PP13-F 

0( 14)=PP14-P 

CF(U=CCNF1 

CF(21=CCNF2 

CF  (3)=CCNF3 

CF(4J=CCNF4 

CF(5)=CCNF5 

CF(6)=CCNF6 

CF  (7)=CCNF7 

CF( 8)=CCNF8 

DC  4911  1=1,8 

CC=CF( I ) 

C(  I,J»  = C( I, JJ+CC 
IF(NWRITE.EC.O)GOTC  4911 
V>PITE(6 ,10) 

wRITE(6,12) I,P, (C (I ,L) ,L=1,9) 
CCNTINJE 

DC  4997  1=1,  NP 

0C=0{ I ) 

U ( I , JJ  =U(  I, J)+0D 

^^  ( I ,J)  =V»(  I , J)+ABS  (CD) 

S(I,J)=S( 1,J)+DD*CC 
IF  (NWRITE.EQ.O)  GCTO  4997 
VvRITE  (6,10) 

WRITE  (6,12)  I,P  , (U( I ,L) ,L  = 1,9) 
WRITE  (6,12)  I,P  , (W{I,L)  ,L=1,9) 
WRITE  (6,12)  I,P  ,(S(I,L)  ,L  = 1,9) 
CCNTINJE 
CCNTINJE 


C 

C 

C 


4999  CONTINUE 


PRINT  SUMMARY  STATISTICS 

DC  5100  J=1 ,9 
51CC  D(J)=.1*J 

WRITE  (6,7)  (0( J ) ,J=1,9 ) 

WRITE  (6,10 
DC  5300  J=l,9 
IF  (NN(J).GT.O)  GCTO  5150 
DC  5125  1=1, NP 
S( 1,J)=1E9 
U(I , J)=1E9 
C( I,J)=1E9 
W( I,J)=1E9 
GCTC  5300 
X J=l. /NN( J) 

DC  5250  1=1, NP 
S(I,J)=SCRT(S(I,J)WXJ) 

U( I , J )=U( I , J) *XJ 
C( I,J)=C( I, J)*XJ 
W(I ,J)=W( 1 , J)*XJ 
CONTINUE 

WRITE  (6,9)  (I,  (U(I,J),J  = 1,9),I  = 1,NP) 
WRITE  ( 6,10 

WRITE  (6,14)(I,(W(I,J),J=1,9) ,I=1,NP) 
WRITE  (6,10) 

WRITE  ( 6,16)(I,(S(I,J),J=1,9) ,I=1,NP) 
WRITE  (6,10) 

WRITE(6,20) (I, (C( I ,J) jJ=l,9 ) ,1=1,8 ) 
FCRMATi  IX,'  C*  ,11,  9F8.3,3X,'  C(;NF'  ) 
wRITE(6,10) 

WRITE  ( 6,15)  (NN( J) ,J=1  ,9) 

IF  (NENC.NE.O)  GOTC  108 
WRITE  (6,11) 

GCTO  25 
ENC 


5125 

515C 

52CC 


52  5C 
53CC 


20 
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